Large scale envrionmental changes may drive changes towards different vegetation types. Biomes with alternative stable states that could be tipped into different vegetation types with small changes in environment create uncertainty in the future vegetation response, and fuel fears of forest instability (Staver et al. 2011). Globally, savannas and forests are often cited to be alternative stable states. However, there is widespread debate over whether alternative savanna and forest states exist or are an artifact of satellite processing (CITE Modis paper, STAVER response), and high uncertainty in predicting their sensitivity to environmental changes (Higgins and Scheiter 2012). Therefore, deterimining if savanna and forests existed as alternative stable states in the recent past and evaulating their past responses to environmental changes can provide key insights into the impact of future changes on ecosystems with alternative states. To explore the potential responses of alternative stable states to environmental changes, we use historical (pre-European Agricultural Settlement) and modern vegetation data from the North American Midwest to ask whether there is evidence for savanna and forest alternative stable states in the temperate zone and determine if these conditions have been resilient to large scale fire suppression, changes in climate.
In the North American temperate zone, much of the savanna-forest boundary has transitioned into a closed forest state (non-agricultural landscape) (Nowaki and Abrams 208). This shift has occured over a time period of large scale fire supppression at least since the 1930's (Nowacki & Abrams 2008). We estimate pre-European agricultural settlement stem density across the Upper Midwest using historical survey data collected in the 19th century (Public Land Survey (PLS)). We evaluate whether there is evidence of alternative stable states in the past by determining if the past landscape had a bimodal distribution of tree density that cannot be accounted for using climate or soils. We then assess how persistant the hypothesized alternative stable states are on the landscape by comparing the historic relationships between density, soils, and climate to the modern relationship of tree denisty to soils and climate. Finally, we apply the past and modern relationships with soils and climate to assess the resilience to projected future climatic changes. We hypothesized that PLS distribution of tree density was bimodal and not well explained by environmental factors due to latent intact vegetation-disturbance feedbacks along the savanna-forest boundary before settlement, but fragmentation and fire suppression has stabilized modern vegetation and produced a modern boundary that is not bimodal.
Understanding forest response to environmental changes is a key challenge to predicting and planning for the future. Vegetation response may be a deteriministic, smooth response to climate in some places, but non-linear and non-deterministic vegetation-climate relationships can result in divergent forest responses. Globally, transistions from savanna to forests are sharp, and often not deterimined soley by climate. Rather, savannas and forests may be “alternative stable states” that coexist within the same climate due to vegetation structure-disturbance feedbacks (Staver et al.2011, CITE others). Savanna and forest alternative states are largely observed in the tropics, fueling fears that the biomes are potentially unstable and that environmental changes and increased disturbances could result in a a collapse of closed forests to the low tree density savanna state (Staver et al. CITE others). The potential collapse of forests that can switch to their alternative savanna stable state would have enormous consequences for terrestrial carbon storage and biodiversity. Similar vegetation structures (savannas and forests) are also observed in the temperate zone, for example, historically, in the North American midwest, the Humid Pampas region in South America, and (temperate steppe in Asia), suggesting a potential that vegetation structure-disturbance feedbacks could have lead to alternative stable states in these regions as well. However, large scale environmental changes (land use, fire suppression, climate changes) have already occurred in the temperate zones, making these regions ideal to study the resilience potential savanna/forest instability to past environmental changes. [paragraph purpose: savannas and forests as alt. stable states, temperate tropics potentially as well, ??uncertainty in envtl. response, we could study this in the temperate zone]
While the climates of temperate and tropical forest/savanna systems differ, the proposed mechanisms for savanna and forest vegetation in both regions converge on disturbance regimes (however, this is well studied in the tropics, where disturbance regimes are intact). While climate does not deterministically predict the vegetation state overall in these regions, precipitation amount, seasonality, and soil factors can all play a role (Lehmann et al 2011, Staver et al. 2011, Danz et al. 2011). Particularly at a large regional scale, deterministic high and low precipitation creates forests and open savs./grasslands. respectively (Danz et al. 2011, Staver et al. 2011). Additionally, hetergeneity in soil factors that affect hydrology and overall climate may be important at smaller spatial scales (Danz 2011). Soil and topography are hypothesized to act as fire breaks, creating upland forested regions and lowland, flat prairies in the temperate zone (grimm, nowacki & abrams). Vegetation structure-disturbance feedbacks are important intermediate precipitation levels in seasonally dry tropical regions (Staver et al. 2011). In these intemediate climates, fire is associated with open savanna ecosystems, where high grassy fuel loads in open systems promote fire spread, which prevents understory trees from establishing. However at these same climates, closed forests may be present but light limition in the understory prevents fuel buildup, suppressing the likelihood of fire (@hoffmann_2012). These vegetation-disturbance feedbacks drive will drive regions of intermediate tree cover either towards a closed forest state in the absence of disturbance, or towards an open savanna state in the presence of disturbance. Thus, savanna and forest stable states in these regions often results in a hallmark bimodal distribution in tree cover that is not deterministically predicted by climate or environmental factors. The instability and stochastic nature of alternative stable states and disturbance feedbacks make it difficult to predict the outcome of these vegetation processes into the future. [Paragraph Purpose: Alternative stable state mechanism and background, temperate and tropical hyp. mechanisms are extremely similar]
In this study, we use historical vegetation data from the North American Midwest to ask whether there is evidence for savanna and forest alternative stable states in the temperate zone and determine if these conditions have been resilient to large scale fire suppression, changes in climate.
*Public Land Survey Data:* The Public Land Survey was commissioned by the U.S. General Land Office to assess the quality of land/timber and assign land titles for Euro-american settlers. The surveyors placed corner posts at 1 mile section corners in a grid within each township. At these section corners, and at 1/4 section corners, the surveyors would record the distance and azimuth from the corner to the nearest two trees, the species of the nearest two trees, and the diameter at breast height (DBH) of those trees. After these historic records were digitized (cite Mladenoff et al. ), we estimate the denisty and basal area at each corner point using the unbiased morista density estimator with correction factors for surveyor bias (@goring_2016, CITE Cogbill). Tree density is averaged across 8km x 8km grid cells in the upper Midwest.
*Forest Inventory Analysis Data:*
Add FIA data methods from Sean/Simon/Jack. FIA estimated tree denisty is also aggregated to the same 8km raster grid resolution.
Climate Data
Mean annual temperature and Mean annual Precipitation data for the modern and historic vegetation-climate relationships are from PRISM datasets. Modern mean annual temperature and precipitation fore each 8km grid cell is calculated from the 800m 30-year Normals dataset. Historical climate is estimated as the average annual Temperature and average annual Precipiation over the 1895-1910 period. While this historical climate period does not overlap with the entirety of the European agricultural settlement era, this climatic period is the oldest for which we have reliable estimates of temperature and precipitation in the region. Precipitation seasonality was calculated from Monthly Prism datasets as follows:
SI = sum 1-12(abs(mi - P/12))/P * 100
Temperature is calculated using the coefficient of variation method TSI = sd(m1….m12)/Tavgannual *100
Soil data
Soils data are derived from statewide gSSURGO 10m raster data product (https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/home/?cid=NRCS142P2_053628). The weighted averages from the top 30cm of the soil were calculated for % sand, available water content (awc), and saturated hydraulic conductivity (ksat) soil parameters using the gSSURGO On-Demand ArcGIS toolbox (https://github.com/ncss-tech/ssurgoOnDemand.git). Weighted averages were calculated on the statewide soil raster datasets for Minnesota, Michigan, Wisconsin, Indiana, and Illinois, then weighted average rasters were mosaicked together to produce one single raster dataset. Maps of soil parameters were then aggregated to a 1km grid scale and an 8km Great Lakes St. Lawrence projected grid scale (PalEON dataset scale).
Pricipal Component Analysis
Since many climate and soils variables are often correlated, a pricipal components analysis was conducted to reduce dimensions and co-linearity in the environmental data. PCA was coducted on the scaled variables of historic MAP, historic Mean temperature, historic Precipitation Seasonality Index, Historic Tempearture Seasonality, AWC, and %sand. The modern principal component scores for climatic and soil variables were predicted using the historic principal component model. This method maintains the PCA of modern and historic environmental data on the same scale. The 1st Principal Component (PC1) was then used as a covariate when assessing bimodality.
Bimodality analysis
Criteria for bimodality uses a both a bimodality coefficient (BC) and Hartigans Diptest for unimodality. Bimodality Coefficient is calculated from size, skewness, and kurtosis of the distributions (Pfister et al. 2013). A BC greater than 0.556 suggests bimodality. Bimodality coefficient (BC) is calculated from the distribuiton of data as follows using the R package ‘modes’:
BC = \(\frac{(skew^2 + 1)}{(kurtosis + 3)*\frac{(n-1)^2}{(n-2)*(n-3)}}\)
Desnity distributions of the PLS and FIA data were estimated. We then calculated the BC and performed Hartigan’s diptest for unimodaltity for the on these smoothed density distribuitons for the PLS and FIA overall, and the BC within different bins of precipitation and the 1st Principal Component. Hartigan’s diptest has a null hypothesis of a unimodal distribution of the data (CITE). This provides one BC value for a given range of precipitation/PC and whether the distribution is significantly different from a unimodal distribution. These ranges are used to determine the climate space where savannas and closed forests coexist.
Generalized Additive Models
In addition to evaluating the BC for ranges of data, we developed several generalized additive models (GAMs) to determine if tree density can be deterministically predicted by climate and/or soils. GAMs were developed with both strictly additive terms and with flexible smooth terms. The full PLS data and FIA data were separated into two random halves that make up the testing and training datasets. GAM’s were developed with the training dataset, and prediction error was estimated using the test dataset. Model parsimony was evaluated basing on AIC values, and model fit was assessed using RMSE.
Future Forest Stability under CMIP5 projected climate scenarios
To investigate the consequences of the modern climate-vegetation relationship, and the past climate-vegetation relatiohsip would have for the stability of future midwestern forests, we utilized climate projections from CMIP5 representative concentration pathways (rcps) 2.6, 4.5, and 8.5. We obtained average monthly precipitation, total annual precipitation, average monthly temperature, and average annual temperature projections for 2070-2099 under the three RCP scenarios from the CCESM4 GCM climate downscaled climate projections (from worldclim, cite). We calculated precipitation and temperature seasonality as described above. Principal componenet scores (based on the past climate PCA) were assigned for each grid cell using the projected climate data from each RCP.
We then identified the climate space (within PC1) that was bimodal in the PLS time, and mapped where this climate space is projected to be in 2070-2100. We also identified the climate space where FIA forests were bimodal on the modern landscape, and where this climate space would be in 2070-2100.
Distribution of modern and past data
The distribution of tree density in the past is significantly bimodal across the entire upper midwest, with a prominant low tree density mode at 30 trees/ha and a high tree density mode at 205 trees/ha (Diptest p < 0.05, BC = 0.74306). While Hartigan’s diptest shows shows a significant multimodality to the Modern FIA tree denstity, the two modes both occur at intermediate and high tree denisty (147 trees/ha, and 550 trees/ha). Therefore, we do not see a signifcant bimodality with low and high tree denisty models on the modern landscape. This shift in tree density is well explained by the original vegetation. Specifically, grid cells with historically low tree density have generally increased in tree denisity, while those that were dense closed forests in the past have decreased in tree density (Linear reg., R^2 = XXXX, p < 0.05, Figure 2). In the past, grid cells with low tree denisty were comprised of mostly Oak and Hickory species. However, grid cells on the modern landscape with low species density generally have mixed species composition of mixed or single species hardwoods. Most low density FIA grid cells are not dominated by Oak or Hickory species, representing a shift in species composition that accompanied the shift in tree denisty.
PCA (supplemental material??)
Principal component 1 (PC1) explains 54% of variance and component 2 (PC2) explains 29% of the variance. Past precipitation and Past Temperature both have strong positive loadings on PC1, while temperature seasonality and precipitation seasonality index both have negative loadings fro PC1. Sand has the strongest negative loading and AWC has strongest positive loading on PC2. Plotting density/vegetatation type in Principal component space does not separate vegeatation into savanna, prairie, and forest.
Generalized Additive Models/Relationship between historic vegetation and climate/soil variables
Gam model selection for the model with the lowest AIC value yeilded the model: tree density ~ s(MAP) + s(mean annual temperatue) + s(deltaT) + s(deltaP). While AIC suggests this is the best model (Table 1), it only explains 61% of the Deviance and has high Mean squared prediction error for predicting the test data set (MSPE = 5976.366). Prediction of tree density is overestimated in low tree denisty grid cells and underestimated in high tree denisty grid cells.
The Gam model for the modern landscape explained only 29 % of the deviance and had a high mean squared prediction error for predicting the test data set (MSPE = XXXX). Prediction of tree density is again overestimated in low tree denisty grid cells and underestimated in high tree density grid cells.
Notes: perhaps we should find the places that there is bimodality and run the gams in these places to see if soils or something else might predict tree density within these bimodality regions.
Historic climate relationship
PLS tree density was significantly bimodal between xx - xx scores of PC1, which corresponded to intermediate values of temperature and precipitation (XXXX mm/year and XX - XX degC). Over 34% (212,032 km^2) of the the past landscape was bimodal (figure 4) during the PLS time.
Modern climate relationship Modern tree denisty was not significantly bimodal overall, and bimodality was detected at less than 1% of the landscape. However, if the historic climate vegetation relationship had remained intact, and only climate had shifted, we would have expected 34.5% of the modern landscape to be bistable today (Figure 4).
Tree density in the past was significantly bimodal, and was not accounted for by climate or soil factors, supporting our hypothesis that savannas and forests were alternative stable states in the region. Projecting the past climate-vegeation relationship onto the modern climate would drive shifts in the climate space where bistable tree distributions are possible, but continue to predict large extent of savanna and forest alternative stable states on the modern landscape. However, modern tree density is no longer bimodal on the modern landscape, suggesting an unexpected shift towards stable closed forests in the region. Neither past tree density nor modern tree density has a fully deterministic relationship with climate, but bimodality in past tree density is predominant at intermediate temperature and precipitation space. These largescale shifts from alternative savanna and forest stable states towards stable closed forests indicate that land use changes (fire suppression, logging, agricultural conversion) have had a far greater impact on midwestern woody ecosystems than changes in climate over the last ~150 years.
Land use changes that have led to fire suppression and land conversion in the region have long been hypothesized to drive “mesophication” (Nowacki and Abrams 2008), and increases in woody cover in the region. Additionally, small scale studies indicate that historic geography that would act as fire breaks (rivers and topography) are important explanatory variables for PLS vegetation (Grimm 1984). Modern land management strategies for restoring open savannas and grasslands in the region require extensive cutting and prescribed burns to limit understory tree growth (CITE). (CITE) suggest that a shift in overall fire return intervals have shifted southern midwest savanna and grasslands towards mesophytic forests and require larger and longer efforts of prescribed burns to return these ecosystems to their open states. While there has been small scale evidence that fire and disturbances were historically important for the maintence of temperate savanna-forest boundary in the US, our study is the first to document evidence for alternative savanna and forest stable states across a large region of the midwest. Additionally, the modern landscape after long term land use changes and post-settlement fire suppression has become mostly stable forests.
Evidence for alternative savanna and forest stable states in the temperate Midwest provides evidence that bistable systems are not constricted to the tropics. In the tropics, intermediate precipitation and high precipiation seasonality is the climatic domain for alternative stable states. However, our finding of temperate alternative stable states at much lower precipitaiton levels and temperatures, indicate that the constraints on growth identified in the tropics do not impose hard limits on the extent of these systems globally. Rather, we propose that alternative stable states are constrained to regions of intemediate moisture, where disturbance feedbacks occur.
Overall the shifts between the modern landscape and the pre-European settlement is marked by large scale forest stabilization, a shift from bimodal distribution in tree cover characterized by closed forests and open savannas, to mostly closed forested ecosystems. Throughtout the holocene, shifts in tree composition and tree cover have been linked to climatic variability across the region (CITE). However, this recent shift is not predicted by the observed changes in temperature, precipitation, or T/P seasonality. Rather, we propose that increased fire suppression and removal of disturbances from the landscape as likely mechanisms for overall shift towards closed forests. However, lack of direct evidence for the mechanism of these shifts makes it difficult to direct attribute the cause. Additonally, increases in atmospheric CO2 have been hypothesized to have a CO2 fertilization effect on tree growth, and could increase forest cover by increasing growth and tree recruitment (Wycoff and Bowers, CITE recruitment paper). To determine how much a CO2 fertilization effect may have contributed to increased forest cover in the Midwest, mechanistic evidence for CO2 fertilization in the region is needed.
The distribution of tree density across the region and the stabiliity of ecosystems to envrionmental changes can have large consequences for future biodiversity, carbon storage, and management decisions. Although evidence for alternative savanna and forest stable states in the midwest and throughout the globe suggest the potential for large scale forest/biome instability with changes in climate, ongoing land use changes and fires suppression have historically stabilized forests in the midwest. While increased stability of terrestrial ecosytems is a positive effects of the “mesophication” and overall conversion of savannas and grasslands to forests, it may come at the large cost of species biodiversity loss, increased vegetation homogenity, and loss of habitat. Assuming current land cover strategies continue, the modern forests are not likely to revert back to savanna in response of climate. However, increased likelihood of fire disturbances and further land use changes that may accompany changes in climate could alter the stability of the region. Being able to identify forest stability, and vegetation responses to different types of environmental changes will be important to understanding the trajectory of future forest and savanna biomes.
Using the climate space with reduced dimensionality (PC1) and the vegetation classificaitons from rheumtella, we can map out places where tree density was bimodal in the past:
Additionally, the biplot for the Principal Components Analysis:
Table 1: PLS density GAM models for stable sites
##
## Family: gaussian
## Link function: identity
##
## Formula:
## PLSdensity ~ s(pasttmean) + s(MAP1910)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 191.284 2.843 67.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(pasttmean) 8.433 8.903 17.55 <2e-16 ***
## s(MAP1910) 7.169 8.237 33.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.344 Deviance explained = 35.2%
## GCV = 10094 Scale est. = 9958 n = 1232
| model | formula | df | AIC | MSPE | dev.expl | |
|---|---|---|---|---|---|---|
| PLS.gam15 | PLS.gam15 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 42.94 | 14186 | 5919 | 63.93 |
| PLS.gam14 | PLS.gam14 | s(sandpct) + s(awc) | 39.43 | 14203 | 5898 | 63.21 |
| PLS.gam13 | PLS.gam13 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 26.47 | 14365 | 6618 | 57.16 |
| PLS.gam12 | PLS.gam12 | s(sandpct) + s(awc) | 24.05 | 14435 | 7124 | 54.48 |
| PLS.gam14L | PLS.gam14L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 6 | 14616 | 8111 | 45.7 |
| PLS.gam15L | PLS.gam15L | s(sandpct) + s(awc) | 6 | 14616 | 8111 | 45.7 |
| PLS.gam12L | PLS.gam12L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 14622 | 8165 | 45.36 |
| PLS.gam11 | PLS.gam11 | s(sandpct) + s(awc) | 27.93 | 14803 | 9093 | 39.03 |
| PLS.gam10 | PLS.gam10 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 20.08 | 14817 | 9304 | 37.54 |
| PLS.gam9 | PLS.gam9 | s(sandpct) + s(awc) | 25.66 | 14853 | 9542 | 36.25 |
| PLS.gam7 | PLS.gam7 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 17.6 | 14857 | 9668 | 35.22 |
| PLS.gam1 | PLS.gam1 | s(sandpct) + s(awc) | 9.717 | 14997 | 10605 | 26.44 |
| PLS.gam13L | PLS.gam13L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 15005 | 10986 | 25.41 |
| PLS.gam10L | PLS.gam10L | s(sandpct) + s(awc) | 5 | 15017 | 10995 | 24.68 |
| PLS.gam11L | PLS.gam11L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 6 | 15018 | 10967 | 24.77 |
| PLS.gam8 | PLS.gam8 | s(sandpct) + s(awc) | 9.296 | 15054 | 11186 | 22.94 |
| PLS.gam9L | PLS.gam9L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 15080 | 11685 | 20.74 |
| PLS.gam7L | PLS.gam7L | s(sandpct) + s(awc) | 4 | 15098 | 11878 | 19.42 |
| PLS.gam2 | PLS.gam2 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 10.68 | 15101 | 12069 | 20.13 |
| PLS.gam5 | PLS.gam5 | s(sandpct) + s(awc) | 5.492 | 15126 | 12040 | 17.76 |
| PLS.gam3 | PLS.gam3 | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.658 | 15127 | 12161 | 18.25 |
| PLS.gam4 | PLS.gam4 | s(sandpct) + s(awc) | 10.15 | 15158 | 12365 | 16.29 |
| PLS.gam8L | PLS.gam8L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 4 | 15175 | 12246 | 14.24 |
| PLS.gam5L | PLS.gam5L | s(sandpct) + s(awc) | 3 | 15180 | 12384 | 13.76 |
| PLS.gam2L | PLS.gam2L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 15180 | 12832 | 13.73 |
| PLS.gam6 | PLS.gam6 | s(sandpct) + s(awc) | 10.5 | 15237 | 13271 | 10.78 |
| PLS.gam4L | PLS.gam4L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 15308 | 14076 | 4.33 |
| PLS.gam6L | PLS.gam6L | s(sandpct) + s(awc) | 3 | 15314 | 14027 | 3.832 |
| PLS.gam1L | PLS.gam1L | PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 15318 | 14176 | 3.572 |
| PLS.gam3L | PLS.gam3L | s(sandpct) + s(awc) | 3 | 15335 | 14007 | 2.196 |
Table 1A: PLS denisty GAM models with log transformed PLSdensity as as the dependant variable
##
## Family: gaussian
## Link function: identity
##
## Formula:
## log(PLSdensity) ~ s(pasttmean) + s(MAP1910)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.88276 0.02488 196.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(pasttmean) 4.055 5.058 34.79 <2e-16 ***
## s(MAP1910) 8.365 8.878 58.83 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.362 Deviance explained = 36.9%
## GCV = 0.77084 Scale est. = 0.76244 n = 1232
| model | formula | df | AIC | MSPE | |
|---|---|---|---|---|---|
| PLS.log15 | PLS.log15 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 43.08 | 2656 | 0.5262 |
| PLS.log14 | PLS.log14 | s(sandpct) + s(awc) | 38.88 | 2677 | 0.5305 |
| PLS.log13 | PLS.log13 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 27.78 | 2725 | 0.5611 |
| PLS.log12 | PLS.log12 | s(sandpct) + s(awc) | 22.97 | 2884 | 0.6378 |
| PLS.log11 | PLS.log11 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 29.54 | 2996 | 0.7545 |
| PLS.log10 | PLS.log10 | s(sandpct) + s(awc) | 24.9 | 3011 | 0.7629 |
| PLS.log12L | PLS.log12L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 3065 | 0.7532 |
| PLS.log14L | PLS.log14L | s(sandpct) + s(awc) | 6 | 3066 | 0.7541 |
| PLS.log15L | PLS.log15L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 6 | 3066 | 0.7541 |
| PLS.log9 | PLS.log9 | s(sandpct) + s(awc) | 23.71 | 3095 | 0.7891 |
| PLS.log7 | PLS.log7 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 14.42 | 3177 | 0.8383 |
| PLS.log10L | PLS.log10L | s(sandpct) + s(awc) | 5 | 3274 | 0.9129 |
| PLS.log11L | PLS.log11L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 6 | 3276 | 0.9128 |
| PLS.log1 | PLS.log1 | s(sandpct) + s(awc) | 10.12 | 3341 | 0.9737 |
| PLS.log9L | PLS.log9L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 3374 | 0.9777 |
| PLS.log13L | PLS.log13L | s(sandpct) + s(awc) | 5 | 3374 | 0.9901 |
| PLS.log8 | PLS.log8 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 16.11 | 3400 | 1.045 |
| PLS.log7L | PLS.log7L | s(sandpct) + s(awc) | 4 | 3434 | 1.023 |
| PLS.log5 | PLS.log5 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 10.6 | 3447 | 1.09 |
| PLS.log8L | PLS.log8L | s(sandpct) + s(awc) | 4 | 3502 | 1.133 |
| PLS.log3 | PLS.log3 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.806 | 3506 | 1.095 |
| PLS.log5L | PLS.log5L | s(sandpct) + s(awc) | 3 | 3509 | 1.139 |
| PLS.log2 | PLS.log2 | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.994 | 3589 | 1.184 |
| PLS.log4 | PLS.log4 | s(sandpct) + s(awc) | 10.55 | 3603 | 1.215 |
| PLS.log2L | PLS.log2L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 3615 | 1.228 |
| PLS.log6 | PLS.log6 | s(sandpct) + s(awc) | 9.716 | 3638 | 1.229 |
| PLS.log6L | PLS.log6L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 3667 | 1.274 |
| PLS.log3L | PLS.log3L | s(sandpct) + s(awc) | 3 | 3677 | 1.278 |
| PLS.log4L | PLS.log4L | log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 3687 | 1.289 |
| PLS.log1L | PLS.log1L | s(sandpct) + s(awc) | 3 | 3717 | 1.33 |
| dev.expl | |
|---|---|
| PLS.log15 | 60.54 |
| PLS.log14 | 59.61 |
| PLS.log13 | 57.21 |
| PLS.log12 | 50.94 |
| PLS.log11 | 46.84 |
| PLS.log10 | 45.81 |
| PLS.log12L | 41.5 |
| PLS.log14L | 41.54 |
| PLS.log15L | 41.54 |
| PLS.log9 | 41.86 |
| PLS.log7 | 36.89 |
| PLS.log10L | 30.66 |
| PLS.log11L | 30.66 |
| PLS.log1 | 27.44 |
| PLS.log9L | 24.86 |
| PLS.log13L | 24.85 |
| PLS.log8 | 24.6 |
| PLS.log7L | 20.98 |
| PLS.log5 | 20.94 |
| PLS.log8L | 16.48 |
| PLS.log3 | 16.98 |
| PLS.log5L | 15.84 |
| PLS.log2 | 11.24 |
| PLS.log4 | 10.27 |
| PLS.log2L | 8.292 |
| PLS.log6 | 7.589 |
| PLS.log6L | 4.355 |
| PLS.log3L | 3.558 |
| PLS.log4L | 2.798 |
| PLS.log1L | 0.3809 |
Table 2: FIA density GAM models
| model | formula | df | AIC | MSPE | |
|---|---|---|---|---|---|
| FIA.gam15 | FIA.gam15 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 31.27 | 19006 | 27066 |
| FIA.gam11 | FIA.gam11 | s(sandpct) + s(awc) | 21.87 | 19298 | 14484 |
| FIA.gam10 | FIA.gam10 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 14.32 | 19321 | 14229 |
| FIA.gam14 | FIA.gam14 | s(sandpct) + s(awc) | 27.44 | 19331 | 26678 |
| FIA.gam10L | FIA.gam10L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 5 | 19346 | 16207 |
| FIA.gam12 | FIA.gam12 | s(sandpct) + s(awc) | 15.18 | 19347 | 25459 |
| FIA.gam11L | FIA.gam11L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 6 | 19348 | 16240 |
| FIA.gam16 | FIA.gam16 | s(sandpct) + s(awc) | 11.96 | 19414 | 14248 |
| FIA.gam14L | FIA.gam14L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 6 | 19421 | 23095 |
| FIA.gam15L | FIA.gam15L | s(sandpct) + s(awc) | 6 | 19421 | 23095 |
| FIA.gam12L | FIA.gam12L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 5 | 19437 | 25197 |
| FIA.gam16L | FIA.gam16L | s(sandpct) + s(awc) | 4 | 19444 | 16088 |
| FIA.gam3 | FIA.gam3 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 10.07 | 19457 | 31573 |
| FIA.gam13 | FIA.gam13 | s(sandpct) + s(awc) | 22.38 | 19462 | 20892 |
| FIA.gam8 | FIA.gam8 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 12.03 | 19490 | 11629 |
| FIA.gam5 | FIA.gam5 | s(sandpct) + s(awc) | 8.802 | 19534 | 11648 |
| FIA.gam8L | FIA.gam8L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 4 | 19542 | 10886 |
| FIA.gam5L | FIA.gam5L | s(sandpct) + s(awc) | 3 | 19547 | 10772 |
| FIA.gam3L | FIA.gam3L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 3 | 19566 | 29092 |
| FIA.gam9 | FIA.gam9 | s(sandpct) + s(awc) | 20.37 | 19645 | 14473 |
| FIA.gam13L | FIA.gam13L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 5 | 19653 | 14621 |
| FIA.gam9L | FIA.gam9L | s(sandpct) + s(awc) | 5 | 19691 | 15974 |
| FIA.gam7L | FIA.gam7L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 4 | 19703 | 16435 |
| FIA.gam1 | FIA.gam1 | s(sandpct) + s(awc) | 8.143 | 19765 | 13264 |
| FIA.gam1L | FIA.gam1L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 3 | 19805 | 18718 |
| FIA.gam6 | FIA.gam6 | s(sandpct) + s(awc) | 10.29 | 19873 | 12076 |
| FIA.gam4 | FIA.gam4 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 7.758 | 19888 | 12145 |
| FIA.gam4L | FIA.gam4L | s(sandpct) + s(awc) | 3 | 19889 | 12690 |
| FIA.gam2 | FIA.gam2 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 9.375 | 19897 | 11659 |
| FIA.gam6L | FIA.gam6L | s(sandpct) + s(awc) | 3 | 19912 | 11451 |
| FIA.gam2L | FIA.gam2L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 3 | 19913 | 11741 |
| FIA.gam7 | FIA.gam7 | s(sandpct) + s(awc) | 16.35 | NA | 14956 |
Table 2A: FIA density GAM models modeled with log transformation of FIA density
| model | formula | AIC | MSPE | dev.expl | |
|---|---|---|---|---|---|
| FIA.gam15 | FIA.gam15 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19006 | Inf | 33.74 |
| FIA.gam11 | FIA.gam11 | s(sandpct) + s(awc) | 19298 | Inf | 19.15 |
| FIA.gam10 | FIA.gam10 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19321 | 5.034e+288 | 17.15 |
| FIA.gam14 | FIA.gam14 | s(sandpct) + s(awc) | 19331 | Inf | 32.78 |
| FIA.gam10L | FIA.gam10L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19346 | Inf | 14.83 |
| FIA.gam12 | FIA.gam12 | s(sandpct) + s(awc) | 19347 | Inf | 31.06 |
| FIA.gam11L | FIA.gam11L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19348 | Inf | 14.84 |
| FIA.gam16 | FIA.gam16 | s(sandpct) + s(awc) | 19414 | Inf | 11.83 |
| FIA.gam14L | FIA.gam14L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19421 | Inf | 26.94 |
| FIA.gam15L | FIA.gam15L | s(sandpct) + s(awc) | 19421 | Inf | 26.94 |
| FIA.gam12L | FIA.gam12L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19437 | Inf | 26.1 |
| FIA.gam16L | FIA.gam16L | s(sandpct) + s(awc) | 19444 | 8.606e+291 | 9.211 |
| FIA.gam3 | FIA.gam3 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19457 | Inf | 25.65 |
| FIA.gam13 | FIA.gam13 | s(sandpct) + s(awc) | 19462 | Inf | 26.56 |
| FIA.gam8 | FIA.gam8 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19490 | Inf | 7.458 |
| FIA.gam5 | FIA.gam5 | s(sandpct) + s(awc) | 19534 | Inf | 4.385 |
| FIA.gam8L | FIA.gam8L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19542 | 2.339e+273 | 3.293 |
| FIA.gam5L | FIA.gam5L | s(sandpct) + s(awc) | 19547 | 4.656e+233 | 2.859 |
| FIA.gam3L | FIA.gam3L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19566 | Inf | 19.62 |
| FIA.gam9 | FIA.gam9 | s(sandpct) + s(awc) | 19645 | Inf | 17.39 |
| FIA.gam13L | FIA.gam13L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19653 | Inf | 15.36 |
| FIA.gam9L | FIA.gam9L | s(sandpct) + s(awc) | 19691 | Inf | 13.27 |
| FIA.gam7L | FIA.gam7L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19703 | Inf | 12.49 |
| FIA.gam1 | FIA.gam1 | s(sandpct) + s(awc) | 19765 | 1.12e+207 | 9.527 |
| FIA.gam1L | FIA.gam1L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19805 | Inf | 6.572 |
| FIA.gam6 | FIA.gam6 | s(sandpct) + s(awc) | 19873 | 2.046e+286 | 3.412 |
| FIA.gam4 | FIA.gam4 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19888 | 1.059e+178 | 2.159 |
| FIA.gam4L | FIA.gam4L | s(sandpct) + s(awc) | 19889 | 6.106e+189 | 1.543 |
| FIA.gam2 | FIA.gam2 | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19897 | 1.182e+210 | 1.819 |
| FIA.gam6L | FIA.gam6L | s(sandpct) + s(awc) | 19912 | 4.357e+166 | 0.06073 |
| FIA.gam2L | FIA.gam2L | FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + | 19913 | 8.102e+165 | 0.009375 |
| FIA.gam7 | FIA.gam7 | s(sandpct) + s(awc) | NA | 7.22e+296 | 17.1 |
The model with the lowest AIC value is represented by the formula: PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + s(sandpct) + s(awc)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -72.75 107.10 199.40 187.10 267.50 437.90
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.664 78.470 174.600 175.500 261.200 757.600
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 85.01 185.50 284.50 278.60 368.70 495.90 114
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 85.01 185.50 284.50 278.60 368.70 495.90 114
The predicted vs. observed plot shows overestimation of PLS tree density in Low density areas and underestimation of PLS tree density in High density areas. Below we map out the predictions and prediction errors in space. Model fit has improved for the stable portions of the PLS landscape, and does an Okay job predicting PLS tree density. Note that the predictions here are overall much higher than the observations for FIA tree denisty.
# Map out predictions in space:
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
FIA models overpredict over the places that have low tree denisty. I think that if we remove the low denisty places, we might be able to better predict.
The above models are predicting density (continious variable) from continous envrionmental and climate covariates. However, we can also use the savanna/forest density classificaitons as our y variable and predict the probability a grid cell is forest (ecocode = 1, PLS density > 47 trees/ha) and savannna (ecocode = 0, PLSdensity < 47 trees/ha & >0.5 trees/ha) from continous environmental and climate covariates.
This method has the benefit of predicting the probability of forest (or the probability of savanna) given certain environmental conditions. For this analysis, we exclude prairie sites
| model | modeltype | formula | df | AIC | |
|---|---|---|---|---|---|
| PLS.lgr14 | PLS.lgr14 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 50.88 | 3704 |
| PLS.lgr12 | PLS.lgr12 | smooth | s(sandpct) + s(awc) | 27.6 | 3875 |
| PLS.lgr13 | PLS.lgr13 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 27.6 | 3875 |
| PLS.lgr11 | PLS.lgr11 | smooth | s(sandpct) + s(awc) | 35.53 | 4697 |
| PLS.lgr10 | PLS.lgr10 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 27.32 | 4808 |
| PLS.lgr9 | PLS.lgr9 | smooth | s(sandpct) + s(awc) | 27.4 | 4997 |
| PLS.lgr14L | PLS.lgr14L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 5248 |
| PLS.lgr12L | PLS.lgr12L | linear | s(sandpct) + s(awc) | 4 | 5257 |
| PLS.lgr13L | PLS.lgr13L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 4 | 5257 |
| PLS.lgr7 | PLS.lgr7 | smooth | s(sandpct) + s(awc) | 18.85 | 5318 |
| PLS.lgr11L | PLS.lgr11L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 5933 |
| PLS.lgr10L | PLS.lgr10L | linear | s(sandpct) + s(awc) | 4 | 6021 |
| PLS.lgr8 | PLS.lgr8 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 18.01 | 6525 |
| PLS.lgr9L | PLS.lgr9L | linear | s(sandpct) + s(awc) | 4 | 6884 |
| PLS.lgr8L | PLS.lgr8L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 6888 |
| PLS.lgr5 | PLS.lgr5 | smooth | s(sandpct) + s(awc) | 9.597 | 6894 |
| PLS.lgr5L | PLS.lgr5L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 6971 |
| PLS.lgr7L | PLS.lgr7L | linear | s(sandpct) + s(awc) | 3 | 7109 |
| PLS.lgr1 | PLS.lgr1 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.147 | 7158 |
| PLS.lgr2 | PLS.lgr2 | smooth | s(sandpct) + s(awc) | 9.813 | 7453 |
| PLS.lgr3 | PLS.lgr3 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.884 | 7645 |
| PLS.lgr6 | PLS.lgr6 | smooth | s(sandpct) + s(awc) | 9.833 | 7823 |
| PLS.lgr4 | PLS.lgr4 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.658 | 7826 |
| PLS.lgr3L | PLS.lgr3L | linear | s(sandpct) + s(awc) | 2 | 7998 |
| PLS.lgr2L | PLS.lgr2L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 8077 |
| PLS.lgr6L | PLS.lgr6L | linear | s(sandpct) + s(awc) | 2 | 8116 |
| PLS.lgr4L | PLS.lgr4L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 8148 |
| PLS.lgr1L | PLS.lgr1L | linear | s(sandpct) + s(awc) | 2 | 8414 |
| MSPE | |
|---|---|
| PLS.lgr14 | 0.09428 |
| PLS.lgr12 | 0.09668 |
| PLS.lgr13 | 0.09668 |
| PLS.lgr11 | 0.1199 |
| PLS.lgr10 | 0.1203 |
| PLS.lgr9 | 0.1246 |
| PLS.lgr14L | 0.1346 |
| PLS.lgr12L | 0.1352 |
| PLS.lgr13L | 0.1352 |
| PLS.lgr7 | 0.1333 |
| PLS.lgr11L | 0.1505 |
| PLS.lgr10L | 0.1497 |
| PLS.lgr8 | 0.1709 |
| PLS.lgr9L | 0.1731 |
| PLS.lgr8L | 0.18 |
| PLS.lgr5 | 0.1792 |
| PLS.lgr5L | 0.1802 |
| PLS.lgr7L | 0.1831 |
| PLS.lgr1 | 0.1762 |
| PLS.lgr2 | 0.1915 |
| PLS.lgr3 | 0.1989 |
| PLS.lgr6 | 0.1993 |
| PLS.lgr4 | 0.2018 |
| PLS.lgr3L | 0.2092 |
| PLS.lgr2L | 0.2095 |
| PLS.lgr6L | 0.208 |
| PLS.lgr4L | 0.2108 |
| PLS.lgr1L | 0.2203 |
Looks like model number 14 with just Mean Annual Precipitation has lowest AIC for the logistic/binomial model. Lets map out the probability of forest based on model with the lowest AIC value:
forest type ~ s(MAP1910)+ s(pasttmean) + s(deltaT) + s(pastdeltaP)
The smooth predictors in this model are the same as the smooth predictors in the denisty model
| model | modeltype | formula | df | AIC | MSPE | |
|---|---|---|---|---|---|---|
| FIA.lgr14 | FIA.lgr14 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 34.44 | 685.7 | 0.1092 |
| FIA.lgr12 | FIA.lgr12 | smooth | s(sandpct) + s(awc) | 22.58 | 717.7 | 0.1066 |
| FIA.lgr13 | FIA.lgr13 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 22.58 | 717.7 | 0.1066 |
| FIA.lgr11 | FIA.lgr11 | smooth | s(sandpct) + s(awc) | 22.15 | 824.1 | 0.1366 |
| FIA.lgr10 | FIA.lgr10 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 20.14 | 832.3 | 0.1367 |
| FIA.lgr9 | FIA.lgr9 | smooth | s(sandpct) + s(awc) | 22.53 | 862 | 0.1458 |
| FIA.lgr7 | FIA.lgr7 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 15.17 | 888.6 | 0.1556 |
| FIA.lgr14L | FIA.lgr14L | linear | s(sandpct) + s(awc) | 5 | 1016 | 0.1397 |
| FIA.lgr11L | FIA.lgr11L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 5 | 1032 | 0.1537 |
| FIA.lgr12L | FIA.lgr12L | linear | s(sandpct) + s(awc) | 4 | 1035 | 0.1406 |
| FIA.lgr13L | FIA.lgr13L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 4 | 1035 | 0.1406 |
| FIA.lgr10L | FIA.lgr10L | linear | s(sandpct) + s(awc) | 4 | 1050 | 0.1514 |
| FIA.lgr8 | FIA.lgr8 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 15.17 | 1085 | 0.1993 |
| FIA.lgr2 | FIA.lgr2 | smooth | s(sandpct) + s(awc) | 9.557 | 1088 | 0.216 |
| FIA.lgr9L | FIA.lgr9L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 4 | 1186 | 0.1772 |
| FIA.lgr4 | FIA.lgr4 | smooth | s(sandpct) + s(awc) | 9.565 | 1192 | 0.2224 |
| FIA.lgr8L | FIA.lgr8L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 1194 | 0.1881 |
| FIA.lgr5 | FIA.lgr5 | smooth | s(sandpct) + s(awc) | 9.793 | 1196 | 0.192 |
| FIA.lgr5L | FIA.lgr5L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 1204 | 0.188 |
| FIA.lgr1 | FIA.lgr1 | smooth | s(sandpct) + s(awc) | 8.648 | 1209 | 0.193 |
| FIA.lgr7L | FIA.lgr7L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 3 | 1234 | 0.1929 |
| FIA.lgr6 | FIA.lgr6 | smooth | s(sandpct) + s(awc) | 9.611 | 1274 | 0.229 |
| FIA.lgr3 | FIA.lgr3 | smooth | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 9.157 | 1297 | 0.2314 |
| FIA.lgr2L | FIA.lgr2L | linear | s(sandpct) + s(awc) | 2 | 1323 | 0.2239 |
| FIA.lgr4L | FIA.lgr4L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 1329 | 0.2207 |
| FIA.lgr6L | FIA.lgr6L | linear | s(sandpct) + s(awc) | 2 | 1400 | 0.2203 |
| FIA.lgr1L | FIA.lgr1L | linear | ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + | 2 | 1434 | 0.2376 |
| FIA.lgr3L | FIA.lgr3L | linear | s(sandpct) + s(awc) | 2 | 1456 | 0.2319 |
## png
## 2
The above map shows ‘savana or forest’ places where the binomial familiy model predicted between 30-70% probability of forest. There are some places where the model predicts deterministic forest and deterministic savanna.
The main effect on predicted forest here is temperature and temperature seasonality.
The past and modern climate space of the Midwest is broad spans a large range of precipitation and temperature.
PLS & FIA tree density are not readily explained by the gradient in Precipitation, but tree denisty bimodality appears to occur between 600-900 mm/year in the past.
PLS & FIA density are not obviously explained as a function of mean annual temperature, but the bimodality in tree density appears to occur between 0-10 DegC in the past:
Ans. % Sand, Available Water Content (awc), and saturated hydraulic conductivity (ksat) do not explain tree density overall on their own. Additionally, PC1 and PC2 from a Principle Components Analysis do not explain tree density.
Ans. We can look at bins of precipitation to get an idea of where this bimodality coefficient is likely to be highest:
We can look at the bins of temperature to get an idea of where the BC is likely to be highest:
We can also look at the density distribution with bins of the first Principal Component of the environmental data:
We can also look at the density distribution with bins of the second Principal Component of the environmental data:
Bimodality coefficient is calculated from a distribuiton of data as follows:
BC = \(\frac{(skew^2 + 1)}{(kurtosis + 3)*\frac{(n-1)^2}{(n-2)*(n-3)}}\)
Here we calculated the BC for a bin range of 25mm of precipitation, with non-overlapping and overlapping bins:
We can also calculate BC of the Principal components (this is not as easily interpretable as binning by precipitation):
Within the places that were bimodal in the past, and that are bimodal now, we want to know which were savannas and which were forests. We use a classificaiton scheme from Rheumtella et al. (), However there are alternate classificaitons that we could use, such as the scheme Chicago Wilderness has used. We could also make our own classification schemes using kmeans clustering.
| Ecotype | Rheumtella | ChicagoWilderness |
|---|---|---|
| Forest | > 47 trees/ha | > 100 trees/ha |
| Savanna | 0.5-47 trees/ha | 10-100 trees/ha |
| Prairie | < 0.5 trees/ha | < 10 trees/ha |
Within precipiation ranges that have bistablity, differences in temperature may account for differences in tree density. Other regions that are bimodal w.r.t. precipitaiton can be explained by soil factors, such as differences in soil sandiness (PC2). However, a small subset of areas are bimodal in both PC2 and PC1 space.